SU(iV) magnetism in chains of ultracold alkaline-earth-metal atoms: Mott transitions 

and quantum correlations 



Salvatore R. Manmana, 1 Kaden R. A. Hazzard, 1 Gang Chen, 1 Adrian E. Feiguin, 2 and Ana Maria Rey 1 

1 JILA, (University of Colorado and NIST), and Department of Physics, 
University of Colorado, Boulder, Colorado 80309-0440, USA 
2 Department of Physics and Astronomy, University of Wyoming, Laramie, WY 82071, USA 

(Dated: October 6, 2011) 

We investigate one dimensional SU(7V) Hubbard chains at zero temperature, which can be em- 
ulated with ultracold alkaline earth atoms, by using the density matrix renormalization group 
(DMRG), Bethe ansatz (BA), and bosonization. We compute experimental observables and use 
the DMRG to benchmark the accuracy ol the Bethe ansatz for N > 2 where the BA is only ap- 
proximate. In the worst case, we find a relative error e < 4% in the BA ground state energy for 
N < 4 at filling 1/A, which is due to the fact that BA improperly treats the triply and higher 
occupied states. Using the DMRG for N < 4 and the BA for large N, we determine the regimes of 
validity of strong- and weak-coupling perturbation theory for all values of N and in particular, the 
parameter range in which the system is well described by a SU(iV) Heisenberg model at filling 1/N. 
We find this depends only weakly on N. We investigate the Berezinskii-Kosterlitz-Thouless phase 
transition from a Luttinger liquid to a Mott-insulator by computing the fidelity susceptibility and 
the Luttinger parameter K p at 1/N filling. The numerical findings give strong evidence that the 
fidelity susceptibility develops a minimum at a critical interaction strength which is found to occur 
at a finite positive value for N > 2. 

PACS numbers: 67.85.-d, 37.10.Jk, 71.10.Fd, 03.75.Ss 



I. INTRODUCTION 

The SU(A r ) Hubbard model describes A^-fLavor 
fermions hopping on a lattice with flavor-independent on- 
site interactions. The model is a generalization of the 
conventional SU(2) Hubbard model introduced in the 
1960s for the theoretical description of itinerant ferro- 
magnets in spin 1/2 systems [1-3]. It has attracted con- 
siderable theoretical attention in recent years [4-33]. In 
the strong interacting regime and at 1 /N filling (one par- 
ticle per site), the low energy physics of the S\J(N) Hub- 
bard model is captured by an effective SU(A^) Heisen- 
berg model in which the charge degrees of freedom is 
frozen and only the spin degrees of freedom are al- 
lowed to fluctuate. There is a long history of studies of 
SXJ(N) spin systems [4, 34-37]. Their initial motivation 
was to better understand the usual SU(2) antiferromag- 
nets since S\J(N) spins are analytically tractable in the 
large- N limit. These studies have found rich phase dia- 
grams [34, 38-52], exhibiting antiferromagnetic ordering, 
valence-bond solids and in ID exotic spin-nematic phases 
[53-61] and generalizations of the AKLT state [37, 62- 
70], among others. However, as no exact SXJ(N) models 
have existed in nature, these predictions were considered 
as a theoretical playground. 

The recent discovery that the SU(A^) Hubbard model 
describes ultracold gases of alkaline earth atoms on opti- 
cal lattices [71-74] brings these considerations into a new 
perspective and has spurred theoretical and experimen- 
tal interest. Having in mind this specific experimental 
implementation, exotic new phases - such as chiral spin 
liquids [75, 76] - have been predicted [67] which can be 
of relevance for the realization of topological quantum 



computers [77, 78]. 

A first step towards the experimental observation of 
the rich spin physics in ultracold alkaline earth atoms is 
the precise knowledge of the parameter regime in which 
the SU(Af) Heisenberg model describes the low energy 
physics of the full SU(A r ) Hubbard model. In this pa- 
per, we address this question. By combining numeri- 
cal density matrix renormalization group (DMRG) [79- 
82] simulations with Bethe ansatz (BA), we are able to 
predict the onset of the validity of the spin-models for 
all values of N. We pursue three main purposes with 
this paper. Firstly, we address the mainly theoretical as- 
pects of the validity of the approximate BA for A^ > 2. 
Secondly, using quantum information measures we deter- 
mine the value of the critical interaction U c at which the 
Mott transition [83] takes place. Thirdly, we use these 
insights to provide predictions of the minimal value of 
U s at which the atoms behave as spin systems. Those 
predictions are relevant for ongoing experiments with ul- 
tracold alkaline earth atoms. We also discuss numerical 
results for the on-site occupancies, the various correlation 
functions, momentum distributions, and their structure 
factors which can be accessed in current experiments. 

We present these results as follows. In Sec. II, we intro- 
duce the S\J{N) Hubbard model and the S\J{N) Heisen- 
berg model and discuss their realization in ultracold al- 
kaline earth atoms on optical lattices. We also describe 
the methods we use to study these models, and their 
connection to experimentally relevant observables. We 
discuss in Sec. II C the BA treatment of the SU( N) Hub- 
bard chains for N > 2, and in Sec. II F the properties of 
the fidelity susceptibility x(U) in the vicinity of a quan- 
tum critical point. In Sec. Ill, we compare the BA and 
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DMRG results for the ground state energy per site in the 
thermodynamic limit and discuss the relative error of the 
BA results as a function of U. In Sec. IV we address the 
Mott transition using our results for the fidelity suscep- 
tibility. Based on the results of these sections, we turn 
to experimentally relevant aspects in Sec. V. In particu- 
lar, we identify U s > I2t for all values of N by compar- 
ing energies and spin correlation functions of the Hub- 
bard and Heisenberg systems in Sees. VA and VC. In 
Sees. VC and VD, we present our DMRG results for the 
on-site occupancies, and the correlation functions and 
their structure factors and discuss possibilities to access 
these quantities in experiments. Section VI summarizes 
our findings. 



II. MODELS, EXPERIMENTAL REALIZATION, 
METHODS AND OBSERVABLES 

A. Models 

We treat the ID fermionic SU(AT) Hubbard model de- 
scribed by the Hamiltonian 



(ij) , a z,a^a' 

(i) 

using Bethe ansatz and DMRG. The fermionic opera- 
tor fa.i (fa i) destroys (creates) a particle of flavor a on 
lattice site i, with a — l,...,N, and Y)u j\ denotes a 
sum over nearest neighbor sites i and j. In the limit of 
large U/t, the effective model obtained by second order 
degenerate perturbation theory is the SU(JV) Heisenberg 
model 

9t 2 

H =Jj £ ( 2 ) 

with the spin operators S^(i) = i / /3 i which satisfy the 
SU(iV) algebra S*(j)] = 6^(6^ - S aS Sj) and 

hence form the generators of the SU(iV) symmetry. We 
consider the deep Mott insulator (MI) state, in which 
the system satisfies the local constraint of having one 
particle per site, J2 a faifai = 1> so that the spin on 
each site transforms in the fundamental representation of 
SU(iV). Thus the number of sites needed to form a singlet 
is N, leading to rich and exotic behavior [37, 84, 85]. In 
Ref. 67 it was shown that on a square lattice, besides the 
expected antiferromagnetic phase, there are valence bond 
solids, and - most interestingly - topologically ordered 
chiral spin liquids. 



B. Experimental Realization in Systems of 
Alkaline Earth Atoms 

Alkaline earth atoms belong to the second column of 
the periodic table and have two-outer valence electrons. 
These and other atoms with similar atomic structure such 
as Yb have unique atomic properties that make them at- 
tractive candidates for new types of atomic clocks [86- 
88], quantum simulation [71, 73] and quantum informa- 
tion applications [89-91]. In their ground state, 1 Sq, the 
electronic degrees of freedom have neither spin nor or- 
bital angular momentum (J = 0) and the nuclear spin 
(I) is thus decoupled from the electronic state. Only 
fermionic isotopes have I > [92] and these are the fo- 
cus of our study. This decoupling not only allows one 
to independently manipulate nuclear and electronic de- 
grees of freedom, but also implies that the 1 S s-wave 
scattering lengths are independent of the nuclear spin. 
The nuclear-spin-dependent variation of the scattering 
lengths is expected to be smaller than ~ 10~ 9 [71]. This 
leads to the SU(iV) symmetric models treated in this pa- 
per, with N < 21+ 1. 

A fundamental consequence of the SU(iV) symmetry is 
the conservation of the total number of atoms with nu- 
clear spin projection a. This means that an atom with 
large I, e.g. 87 Sr (7 = 9/2), can reproduce the dynamics 
of atoms with lower / if one takes an initial state with 
no population in the extra levels. This feature of SU(iV) 
symmetry is in stark contrast to the case of weaker SU(2) 
symmetry exhibited in alkali atoms, where the depen- 
dence of scattering lengths on the total spin of the two 
colliding particles allows for spin changing collisions, due 
to the finite hyperfine interactions. 

The many-body Hamiltonian that describes cold 
fermionic alkaline-earth atoms in the So state loaded 
in the lowest band of an optical lattice is Eq. (1) with 
t = — J d 3 rw(r)(-|gV 2 + V(r))w(r - r ), where r 
is the separation of two nearest neighbors, and U — 
^j-agg f d 3 rw 4 (r), with M the mass, a gg the nuclear 
spin independent 1 S scattering length and w(r) the 
Wannier functions of an atom in a lattice potential V(r) 
[71] . In order to allow motion only along one direction the 
lattice confinement along the other two must be strong 
to suppress tunneling in the course of the experiment. 



C. Approximate Bethe Ansatz for N > 2 

The Bethe ansatz exactly solves many models, for ex- 
ample the SU(2) Heisenberg model, SU(iV) continuum 
fermions, and the SU(iV) Heisenberg model [34, 93]. 
However, the SU(iV) Hubbard model is not solvable with 
standard Bethe ansatz techniques: since lattice sites may 
be occupied by more than two particles, it is impossible 
to reduce the many-particle problem to two-particle scat- 
tering events. Haldane and Choy developed a "natural" 
generalization of the SU(2) Hubbard model Bethe ansatz 
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equations to SU(iV) symmetry, given in Refs. 94 and 
95. The obtained solution, however only approximates 
the Hubbard model, as discussed below. Appendix A 
presents the details of the calculation. 

As shown there, we obtain for the ground state energy 
per site 



fco 

E BA = — 2t I dk cos(k)p c (k), 

-fcn 



(3) 



with the pseudomomenta k and charge rapidity distribu- 
tion p c {k)i where ko is determined by the density via 



the approximate Bethe ansatz is nevertheless exact due 
to the U 2 /t prefactor. 

Additionally the Bethe ansatz is exact when no con- 
figurations have site occupancies with greater than two 
particles per site. This includes the 1/N filling U = oo 
(hard core) limit where it reproduces the exact BA so- 
lution of the SU(iV) Heisenberg model, and the dilute 
limit (n) <C 1 where it reproduces the behavior of SU(A) 
(5-function interacting particles. 



D. Details of the DMRG calculation 



dk p c (k). 



In this paper we have restricted to the balanced case 
with equal population of each of the N spin components, 
where the Bethe ansatz equations simplify to [96] 



p c {k) = hcos(fc) / dk' p c (k')G N (sm(k) - sin(fc')), 

dw c -^ c -^H/4 sinh W -}) (4 ) 



G N (x) 



1 

2^ 



sinh (NUw/A) 



Even though we work with a balanced gas, for our pur- 
poses it is equally convenient to simply solve the equa- 
tions Eq. (A3) numerically, as described in Appendix A. 

While there have been numerous theoretical studies of 
this approximate Bethe ansatz [94-100], no precise quan- 
tification of its accuracy was available. In Sec. Ill we 
provide such a quantification by comparing to numerical 
DMRG results. 

Haldane and Choy have argued that the Schrodinger 
equation is exactly satisfied only for configurations where 
less than three particles occupy a site [97]. This is the 
reason the BA is approximate. 

The simplest quantification of the approximate nature 
of the SU(iV) Bethe ansatz solution is found by consider- 
ing the SU(A) Hubbard model's three particle problem. 
In particular [97] , one finds that 

(x|ff-£(k)|vfc) = ^/(k)(x|P 3 |k), (5) 

with 
/(k) = 



cos 



COS 



and defining the energy E(k) = — 2t[cos(fci) +cos(fc2) 
+ cos(fc3)], the three-particle per site projection opera- 
tor P 3 = S XlX2 5 X2X3 , the state |k) to be the three particle 
Bethe ansatz wavefunction characterized by pseudomo- 
menta fci, ki 1 and k%, and |x) the state with particles at 
positions x\, X2, and X3. It is illustrative here, however, 
to note that even when U — 0, and thus the number of 
configurations with triple or larger occupancies is large, 



Due to the large on-site Hilbert spaces, the DMRG is 
restricted to N < 5. In addition, for the gapless systems 
treated in the following, a large number of density-matrix 
eigenstates m must be kept. We therefore treat systems 
only up to L = 216 lattice sites and keep up to m = 4000 
density matrix eigenstates. We discuss results obtained 
with open boundary conditions (OBC) since the DMRG 
is most efficient in this case. 

Note that an additional restriction appears when com- 
puting the fidelity susceptibility %([/) [Eq. (13)] discussed 
below. The fidelity F(U) [Eq. (12)] is very close to one, 
so that 1 — F(U) ~ 10~ 5 . It is therefore necessary to 
achieve the corresponding convergence in the energy and 
a discarded weight which is significantly smaller than this 
number. These requirements allow to compute xW) reli- 
ably only for systems with up to L = 192 sites for N = 2 
and L = 48 sites for N — 3, keeping up to m = 4000 
and performing 10 sweeps. Additional results for larger 
system sizes and for N = 4 show the same qualitative 
features, but we will not use them for the finite size ex- 
trapolations since the convergence of these calculations 
does not match the requirements for a reliable analy- 
sis of x(U)- Note that the use of non-abelian quantum 
numbers [101-104] might help in future studies to realize 
larger system sizes and treat larger values of N. 



E. Correlation functions and their structure factors 

In this section we introduce the correlation functions 
which we are going to discuss in more detail in Sec. V C. 
The goals are to compare the spin correlation functions of 
the Heisenberg and the Hubbard systems when varying U 
in order to identify the Heisenberg regime of the Hubbard 
chains, and to provide the structure factors of spin and 
charge correlation functions since these are accessible to 
experiments as discussed below. 

For the Heisenberg model, we compute directly 



(6) 



Note that due to the SU(A) symmetry the correlation 
functions along the spin quantization axis and perpendic- 
ular to it are identical, so that it is sufficient to consider 
only Eq. (6). The situation would, however, be different 
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in the presence of symmetry-breaking external fields or 
population imbalance. 

For the Hubbard chains, by expanding the spin opera- 
tors in terms of the Fermi operators f a ,i, and taking 
into account the SU(-/V) symmetry, one finds that the spin 
correlation functions and the associated structure factor 
can be obtained from the difference of density correlation 
functions of two identical and two different flavors of the 
particles, 



S(l,m) = <nf<) - {nfni) 



S{k) 



-^S(l,m) e lk ^- m \ 



(7) 



where n? 



f 1 f 



Note that here summation over 



repeated indices is not implied. 

Complementary to this, density correlation functions 
and their structure factors are 

N a , a (l,m) = <nf<) - (nf) «) 

A/;, Q (fc) = \Y, N «Ah m ) e ifc( '- m) , (8) 

l.m 



and 



N at0 (l,m) = (nfni) - (nf) (ni) 



(9) 



the momentum distribution function n a (k) can be ac- 
cessed via time of flight measurements [111, 112], the spin 
structure factor deep in the Mott insulator phase S(k) is 
accessible via noise correlations in the time of flight mea- 
surements [113], and the density structure factor Af(k) 
can be measured using Bragg scattering [112, 114, 115]. 
By applying a magnetic field, the nuclear spin states can 
be spectroscopically distinguished in transitions to elec- 
tronic excited states. The reason is that the Lande g- 
factor of the excited state significantly differs from that of 
the ground state (e.g. ~ 60% for strontium, Sr [116]) and 
in a biased magnetic field, the various Zeeman transitions 
have different resonant frequencies, as demonstrated in 
Ref. 117. Hence Bragg scattering with light with fre- 
quency near a resonance for state a measures Af a ,a(^)- 
In the balanced SU(iV) gas, from the knowledge of Af(k) 
and J\f a .a(k) one can compute Af a ,/3- in the more general 
- possibly spin imbalanced - case (not considered in the 
present manuscript), one can directly measure N a p. To 
accomplish this, one uses probe light with a frequency ui 
where more than one spin species has response. Tuning 
uj tunes the relative response of the different spin flavors, 
and thus tunes the correlations that are measured by the 
probe, in a well-characterized way. Ref. 118 analyzes the 
SU(2) case in detail. 



Fidelity and Fidelity Susceptibility in the small 
U limit 



for the correlations between particles of the same and 
of two different species, respectively. Again, due to the 
SU(iV) symmetry, it is sufficient to restrict to two ob- 
servables: one with a — ft, and one with a ^ ft (these 
are otherwise independent of a and (3). In addition, it 
is useful to introduce the correlation function and the 
structure factor of the total density, 

N(l,m) = (iV ; total A^ otal ) - (7V; total ) (AC tal ) 

AA(fc) = i^7V(Z,m)e ifc ('- m ), (10) 

l.m 

with jV* otal = J2 a n t- Further information is provided 
by the one-particle density matrix (OPDM) and the mo- 
mentum distribution function, 



ik(l— m) 



(ii) 



Note that due to the SU(iV) symmetry Qa,e(hj) = for 
a 7^ ft so that it is sufficient to consider only g a . a (l, m). 

While for the experiments on optical lattices it is pos- 
sible to measure the correlation functions in real space 
using in situ techniques [105-110], the structure factors 
in momentum space are easier to access. In particular, 



From Bethe ansatz it is well known that for N = 2 
at 1/N filling there is a Berezinskii-Kosterlitz-Thouless 
transition in the charge sector from a gapless metal- 
lic (Luttinger liquid) phase to a gapped Mott-insulating 
phase at U = [83, 119, 120]. While there are strong 
indications that the transition happens at a value U c > 
for N > 2 [5], Ref. 16 questioned this based on quantum 
information measures computed from DMRG, and a sce- 
nario in which U c is either zero or very close to zero was 
proposed. Here, we investigate the behavior of the sys- 
tem at small values of U by computing the fidelity which 
we define as the overlap between ground states at neigh- 
boring points of the coupling constants (here the on-site 
interaction U) [121] 



F{U) = \(MU)\MU + dU))\ 



and the fidelity susceptibility 



X(U) 



2[1--T(EQ] 
LdU 2 



(12) 



(13) 



also known as the fidelity metric [122]. For many phase 
transitions x is expected to diverge in the thermodynamic 
limit (TL), and it has been shown that it possesses a clear 
signature of the critical point already for rather small sys- 
tems [123-125]. However, in Ref. 126 the singular part 
of the fidelity metric in the vicinity of quantum critical 
points was analyzed by a general scaling argument, and 
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it was found that the singular part of the fidelity suscep- 
tibility does not necessarily diverge at a critical point. 
Instead, it can have a minimum at the critical point. 
For the SU(iV) Hubbard model at l/N filling, there is 
spin-charge separation and each sector is described by 
a Luttinger liquid theory. For two independent theories 
the fidelity factorizes and the fidelity susceptibility is ad- 
ditive [125], which leads to the relation 



x(U) = Xp {u) + xAu) 



(14) 



holds. Since the spin sector realizes a Luttinger liquid 
for all values of U, we can safely reproduce the analysis 
of Refs. 125 and 127 which leads to the relation 



1 (d log (K a (U)) 



dU 



(15) 



Here, K a (U) is the spin Luttinger parameter, and due 
to the SU(iV) symmetry, K a (U) = 1 for all values of 
U > and N. In the Luttinger liquid region, Eq. (15) 
holds also for the relation between Xp(U) and K p (U). 
x{U) computed via Eq. (13) reveals the behavior of the 
charge sector and can be applied to investigate the Mott 
transition. As discussed in more detail in Sec. IV, we 
find numerically that x{U) is minimized at the phase 
transition. 



III. COMPARISON BETWEEN BETHE 
ANSATZ AND DMRG 




FIG. 1. (Color online) Bethe ansatz results for the energies 
in the thermodynamic limit for N = 2, 3, 4, 5 (from top to 
bottom) for density n = 1/2 (1/(2JV) filling, top) and n — 1 
(l/N filling, bottom). The dashed lines indicate the weak 
coupling and the strong coupling limits. 



Figure 1 shows the BA results for the ground state en- 
ergy per site in the thermodynamic limit up to N = 10 
for one atom per site n — 1 (l/N filling) and for half 
an atom per site n = 0.5 (1/(2N) filling). Due to the 
expensive numerics, we consider DMRG results only for 
N < 4 and the same values of the filling. Both DMRG 
and the BA show the same qualitative behavior for all 
N: for large values of U, the energy asymptotically ap- 
proaches a constant, while for small values of U it is 
proportional to U. This suggests the presence of two dif- 
ferent regimes and a crossover region or phase transition 
between them. We come back to this point in Sec. V 
where we discuss in more detail the parameter regime 
in which the systems behave as SU(iV) Heisenberg spin 
chains. Note that at l/N filling, n = 1, upon increas- 
ing N, the energies quickly approach an asymptote, so 
that the curves for N = 3 and N = 4 in Fig. 1 are basi- 
cally indistinguishable. This shows the particles become 
effectively distinguishable quickly for density n = 1. As 
discussed in Ref. 96, in the limit N — > oo a generalization 
of the solution of the Lieb-Liniger equation is obtained, 
so that the SU(iV) Hubbard chain in this limit can be 
regarded as a generalized continuum boson system. 

Figure 2 presents the relative difference between the 
DMRG and the BA results for the ground state energies 
per site in the thermodynamic limit for N < 4. We find 
that the relative error for N — 4 is < 4% at l/N filling, 



n = 1, and < 0.7% at l/(2N) filling n = 0.5. For N = 2, 
the relative error is of the order of 10 -4 or smaller and is 
hence not visible on the scale of the plot. This is expected 
since here the BA is exact, and the DMRG is known to 
be capable of obtaining the ground state energy for finite 
systems with a relative error of 10 -6 or better [81, 128]. 
For N > 2, the BA becomes exact for U = and in 
the limit U — > oo, as explained above in Sec. II C. We 
hence expect the absolute errors to be maximal for some 
intermediate value of U. As can be seen in Fig. 2, the 
maximal relative error is obtained at U ~ 3i. The upturn 
at large- U is discussed below. 

The fraction of sites with three or more particles per 
site can be estimated in the non- interacting limit. For 
densities n < 2, interactions suppress the number of 
triple and higher occupancies, so that the non-interacting 
limit yields an upper bound to the number of such config- 
urations, and thus an upper bound of the corresponding 
error. 

In the non- interacting limit, each of the N flavors is 
independently occupied on a site with probability n/N , 
so that the probability of having to particles per site is 
m C N (n/N) m (l - n/N) N ~ m , with m C N = N\/[m\(N - 
to)!]. Thus the probability of having three or more atoms 
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FIG. 2. (Color online) Relative error between DMRG and 
Bethe ansatz for N = 2, 3, 4 (from top to bottom) at densi- 
ties n = 0.5 (1/(2AQ filling, top) and at n = 1 (1/N filling, 
bottom). 




FIG. 3. (Color online) Triple or higher occupancy (P>3) ver- 
sus N for various densities (n — 0.4, . . . , 1 .0, bottom to top 
in steps of 0.2). 



fluctuations to three or more particles per site with am- 
plitude O ((t/U) 2 ), so Eq. (5) suggests an error in the 
energy of the order 0(t). This absolute error becomes 
small in absolute terms as t — > 0, but the relative error 
diverges since the exact energy scales as t 2 /U. Fig. 2 
shows this upturn at large-t/ of the relative error. 



on a site is 



IV. FIDELITY SUSCEPTIBILITY AS A PROBE 
FOR THE MOTT TRANSITION 



P> 3 (N) =1-J2 m C N (n/N) m (l - n/N) N ~ 



N 



= 1 - (1-n/N) 

[n 2 (N -2)(N-1) + 2n(N - 2)N + 2N 2 
K 2(N-n) 2 



.(16) 



This function monotonically increases with N and con- 
verges to its maximum as N — > oo, given explicitly by 



P> 3 (iV 



oo 



1 



-(n 2 + 2n + 2). 



(17) 



Fig. 3 plots the probability of triple and higher occupa- 
tion P> 3 (N) for N — 3 and N = oo as a function of N 
for various densities. We see that at density n = 1, the 
fraction of triple and higher occupancies for N = 3 is 
P> 3 (3) = 0.037 and for N = oo is P> 3 (N = oo) = 0.080. 

Equation (5) suggests that the relative error is on the 
order of the fraction of the triply or higher occupied sites 
times U 2 /(4£), approximating /(k) ~ 1, which is a typ- 
ical value, although for fci = &2 = &3 = 7r/2 it diverges. 
For N = 4 at n = 1, the quantity P> 3 [U/{U)] 2 is .003 
and .03 for U = t and U = 35£, in agreement with the 
results of Fig. 2 (bottom), about .005 and .04. Since this 
estimate grossly overestimates fluctuations in the strong 
coupling limit, we also estimate the fluctuations there 
for n = 1. Second order perturbation theory in t/U gives 



In this section, we discuss the numerical results for the 
fidelity and fidelity susceptibility for systems with N < 4 
at 1/N filling obtained via Eqs. (12) and (13). We start 
by summarizing general considerations and former results 
for the Hubbard chain. 



A. Fidelity susceptibility at U — 0: exact results 

For small values of U, the fidelity susceptibility x(U) 
[Eq. (13)] can be obtained from standard perturbation 
theory. For the SU(iV) Hubbard chain one obtains [122, 
127] 



i,a>0 



(18) 
(19) 



where E n and \n) are the eigenenergies and correspond- 
ing excited eigenstates of the SU(iV) Hubbard chain, and 
-E and |0) are the energy and eigenstate of the ground 
state. 

At the non-interacting point U = 0, the ground state 
susceptibility can be computed exactly. In momentum 
space, we obtain the result in the thermodynamic limit 
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N 


X 


2 


l/(247r ) ~ 0.00422172 


3 


0.0109003 


4 


0.0227492 


5 


0.0416842 


6 


0.0696197 



TABLE I. Numerical values of the fidelity susceptibility \ at 
U = for N = 2, . . . , 6 as obtained from Eq. (20). The 
N = 2 case is also shown in Ref. 129. 



[127], 



x(u = o) 



N(N — 1) 
2(2tt) 3 



f 7T 1 



-7T — 7T — 7T 



& dfc' dg ^(l-n k+q )n k ,(l-n k , 7 ) m 

{ € k+q — e k + e k'-q — 



with the single particle dispersion of non-interacting 
fermions e k = — 2icos(/c) and n k = S(cf — £fc)i with 
6f the Fermi energy. The resulting numerical values of 
X(U = 0) for N < 6 are listed in Tab. I. Note that for 
-/V = 2 a finite value is obtained, indicating that there is 
no divergence of x{U) at the metal-insulator transition. 
In addition, perturbation theory shows the derivative of 
x(U) at U = is negative for N > 2, demonstrating that 
a minimum is to be expected at some fine value of U. As 
we will see next, the critical point in the Hubbard chain 
is indeed characterized by such a minimum of x(?7). 



B. Scaling behavior and nature of \{U = U c ) for 
SU(iV) Hubbard chains at 1/N filling 

Reference 127 analyzed the scaling behavior of \ f° r 
the SU(2) Hubbard model at 1/N filling. As suggested 
by the exact result at U = for N = 2, x{U) is found 
not to diverge at the metal-insulator transition. Apply- 
ing Ref. 126's scaling argument to the SU(iV) Hubbard 
model, we find that the singular part of the fidelity sus- 
ceptibility goes to zero as one approaches the critical 
point: the scaling exponent of the fidelity susceptibil- 
ity near the critical point is given by 2A — 2z — 1 (A = 2 
is the scaling dimension of the Hubbard interaction and 
z = 1 is the dynamical exponent). Therefore, the sin- 
gular part of the fidelity susceptibility vanishes as one 
approaches the critical point. Moreover, a large scale 
Quantum Monte Carlo calculation for the SU(2) model 
in Ref. 129 finds that x(f7) has a local minimum at the 
critical point U c = 0. This indicates that the regular 
part of the fidelity susceptibility behaves rather flat in 
the vicinity of U c = 0. 

The transition for N > 2 at 1/N filling is believed to 
be of the same type as for N = 2, but at a finite value of 
U [5]. Hence, it is natural to expect that also for N > 2 
the transition is identified by a local minimum of x{U)- 
This is further corroborated by the fact that K p and the 
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FIG. 4. (Color online) Finite size extrapolation of \(U = 0) 
for N — 2, 3, 4. The horizontal lines show the exact values of 
Tab. I. For N = 3 and N = 4 we have not taken into account 
the results for the largest system sizes [L = 96 and L — 36, 
respectively) due to a discarded weight of similar order of 
magnitude as 1 — J~{U). The black lines show a quadratic fit 
for N = 2 and N = 3 and a linear fit for N = 4. 



coupling of umklapp terms obey similar RG flow equa- 
tions (of BKT type) for any value of N, so that at 1/N 
filling it follows from Eq. (15) that x should have similar 
behavior for all N near the critical point. Therefore, we 
expect a minimum of x a t U c for all values of N. Note 
that the numerical computation of x{U) is independent 
from the computation of K pi which could show anoma- 
lous behavior at U = [120, 130], and also independent 
from the computation of the charge gap, from which it is 
difficult to obtain accurate values of U c due to the expo- 
nential behavior at the BKT-type transition. 



C. Numerical results 

In order to identify U c , we have computed the fidelity 
F (U) = \{ip (U = 0)\ip o {U))\ for systems L < 48, and 
find no signature of a phase transition in the finite size 
data. We associate this to the BKT nature of the phase 
transition and expect that a discontinuity at U c should 
appear after extrapolating to the thermodynamic limit. 
Due to the large numerical effort associated with such an 
analysis, we refrain from doing so and focus instead on 
the behavior of the fidelity susceptibility. 

We start our discussion by estimating the accuracy 
of our numerical results by comparing to the exact re- 
sults at U = 0. In Fig. 4 we show our finite-size-scaling 
analysis for x{U) and the comparison to the exact re- 
sults of Tab. I. We obtain X (U = 0) N=2 « 0.0043 (ex- 
act value: X ~ 0.00422), X (U = 0) N=3 0.0112 (exact 
value: x » 0.01090) and x(U = 0) N=4 « 0.02670 (exact 
value: x ~ 0.022749). Fig. 4 also shows that for N = 3 
the results for L = 96 lead to a bad extrapolation, and for 



N = 4 the results for L = 36 are not accurate enough for 
our considerations. We therefore restrict the finite size 
scaling to TV = 2 and TV = 3, for which the relative error 
of the numerical results at U = is < 3%, and discuss 
the qualitative behavior of the finite size data for TV = 4. 

Figure 5 shows x(U) f° r systems up to L = 192 (TV = 
2) and L = 48 (TV = 3), and the extrapolation to the 
thermodynamic limit. The finite size results for TV = 4 
show qualitatively similar behavior. Interestingly, the 
results in Fig. 5 show various similarities between TV = 2 
and TV = 3. In particular, in the thermodynamic limit, a 
minimum is obtained at U min — for TV = 2 and U m - U1 « 
X.5t for TV = 3, followed by a maximum at U max ~ 1.6i 
for TV = 2 and J7 max « 3.2* for TV = 3. Note that the 
numerical values of x(U) at the minimum and at the 
maximum for TV = 2 and TV = 3 are very similar to each 
other. This raises the question if these values might be 
universal for all values of TV. 

Finite size results for TV = 4 for L < 36 indicate sim- 
ilar behavior, with U min ~ 2t and U max « 3.5i, but the 
more difficult convergence inhibits obtaining the values 
of x(U) at the minimum and the maximum in the ther- 
modynamic limit. Noteworthy is also the finding of a 
universal crossing point in Fig. 5 between the minimum 
and the maximum of x(U), whose explanation lies be- 
yond the scope of the present paper. We therefore leave 
these issues open for future research. 

The findings of Fig. 5 are interesting. As suggested by 
the analysis of Sec. IV B, they support values of U c = 
for TV = 2 (in agreement with the exact BA result), and 
U c « 1.5i and U c « 2t for TV = 3 and TV = 4, respectively. 
These values of U c can be contrasted to the findings of 
Ref. 5 in which the analysis of numerical QMC results 
for the charge gap indicate U c « 2. It for TV = 3 and 
U c ~ 2.8t for TV = 4. In agreement with the conclusions 
of Ref. 16, this indicates that the analysis of the charge 
gap tends to overestimate the value of U c . However, our 
results for x(U) for TV = 2 and TV = 3 indicate that the 
dependence of U m i n on the system size for this quantity 
is rather weak, so that the extrapolation appears to be 
well controlled. We believe therefore that the analysis 
of Ref. 16 underestimates the values of U c , and that in- 
deed U c > t for TV > 2. This is further corroborated 
by computing analytically xiP) m the limit U — > us- 
ing perturbation theory. For TV > 2, we find that x(U) 
is finite and decreases with U , supporting a scenario in 
which the minimum is located at a finite value of U . Note 
that the values of f7 m i n are in rough agreement with the 
results of an approximate BA treatment of the metal- 
insulator transition [99], which finds 2.5i <U C < 3.5i for 
TV = 3, ... , oo, overestimating the value of U c for TV = 3. 

We complement these considerations by another esti- 
mate of x{U) as obtained from Eq. (15), which relates 
x{U) to K P (U) for the Luttinger model. Numerically 
computing the derivative of K p {U) shows that x(Z7) has 
indeed a minimum which is located at U « for TV = 2, 
U « l.lt for TV = 3 and U « 2. It for TV = 4. These 
values are in good agreement with the values of U m i n ob- 
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FIG. 5. (Color online) Fidelity susceptibility for different sys- 
tem sizes for TV = 2 and TV = 3. The black continuous line 
shows the value after extrapolating to the thermodynamic 
limit using the system sizes shown. 



tained by directly computing x(U), and we associate the 
discrepancy to the errors in the numerical computation 
of the derivative. 

We finish this section by relating our findings to ongo- 
ing experiments. Due to the smallness of the charge gap, 
it will be difficult to precisely locate U c in the experi- 
ments. However, our analysis suggests that for all TV > 2 
Luttinger-liquid (LL) behavior can be addressed by the 
experiments in the regime U < t, and that the transition 
to a Mott-insulator takes place at some value of U < At 
for Af — > oo. This suggests that for TV > 2 it should be 
possible to realize the LL by tuning U to a small enough 
value, but not necessarily exactly to zero. Also note that 
similar considerations imply that it should be possible 
to realize at finite temperatures T spin-incoherent LL 
phases which we expect for hus/L <C ksT <C huc/L, 
with us and uc the velocity of the spin and charge exci- 
tations, respectively [131], which for small (/<t behave 
as us ~ v F - U/(2it) and u c ~ v F + U (TV - l)/(2n), 
where the Fermi velocity vf = 2t sin(7r/TV). 

In the following section we will focus on the strong cou- 
pling limit of 1/TV filled SU(TV) Hubbard chains, which 
maps to SU(TV) spin chains for U large enough. From the 
considerations in this section, we expect this mapping to 
work for U 3> 4i. In the following, we provide a more 



precise estimate by comparing energies and correlation 
functions. 



V. HEISENBERG LIMIT OF SU(7V) HUBBARD 
CHAINS 

A. Energies 

It is a priori unclear for which values of U/t the SU(iV) 
Hubbard model behaves as a Heisenberg model at low en- 
ergies, especially in the light of the probable differences 
in U c for N — 2 and N > 2 discussed in the previous 
section. We address this by first comparing the DMRG 
energies of the SU(7V) Hubbard models to the expected 
~ t 2 /U behavior in Fig. 6 and find that, for the values of 
N shown, the Heisenberg regime starts at Us ~ lit- At 
this value, the difference of the DMRG results for = 2 
to the expected t 2 /U behavior is e w 5 x 10~ 3 . Note that 
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FIG. 7. (Color online) Black lines (top two curves in each 
panel) indicate where relative error of the ground state energy 
is within 1% (lower black line) and 3% (upper black line) of 
the t 2 /U asymptote. Blue lines (bottom two curves in each 
panel) indicate where the relative error of the ground state 
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line) of the weak-coupling perturbation theory. The red curve 
(extra middle line in bottom panel) indicates the critical U c 
for the Mott transition within the Bethe ansatz. Top: 1/(2N) 
filling (n = 0.5), Bottom: 1/N filling (n = 1). 
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FIG. 6. (Color online) DMRG results for the absolute value 
of the ground state energy per site en in the thermodynamic 
limit. Top: weak coupling regime. The black lines are linear 
fits to the energy using the first two data points. Bottom: 
crossover region to the strong coupling regime (log- log scale). 
The black lines are fits to a function ~ t 2 /U using the last 
two data points. 



Us shows a slight decrease upon increasing N indicating 
that for all values of N and for U > lit the system be- 
haves as a Heisenberg chain. This is further confirmed by 
BA, which shows that for N < 10 the energy follows the 
t 2 jU behavior in this regime, as shown in Fig. 7. Note, 
however, that the BA shows a slight increase of Us with 
N. The behavior of the energies hence suggests that in 
ID, for all values of N, SU(7V) Heisenberg physics can be 
realized in the experiments with ultracold alkaline earth 
atoms for U > lit. We will further refine this in the 
next section, where we compare the numerical values of 
spin correlation functions of both models as a function 
of U/t. Note that this numerical value of Us is in good 
agreement with the findings of Ref. 132 for a frustrated 
2D system, in which effective spin models are found to de- 
scribe the SU(2) Hubbard model on the triangular lattice 
for U > lOt. We therefore expect that the SXJ(N) Heisen- 
berg model may quantitatively describe experiments with 
alkaline earth atoms in optical lattices for U > Us also 
in higher dimensions. 
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B. Higher local occupancies 

One specific property of the SU(iV) Hubbard systems 
is that each site can be populated with up to N parti- 
cles. In the context of a possible realization of SU(7V) 
Heisenberg physics, it is therefore interesting to analyze 
the strong-coupling behavior of the higher local occu- 
pancies. These quantities are accessible in experiments 
and can provide valuable information. More specifically, 
measurements of the P m , the number of particles on sites 
with occupancy m, will allow to obtain all the moments 
of the density as well as related quantities of interest such 
as the photoassociation rate [133], which is oc (n(n — 1)) 
(with n — ~^2 a n a ), and the rate at which atoms are lost 
from the trap, dominated by 3-body losses and hence 
oc (n(n — l)(n — 2)). In situ single-site resolved mea- 
surements directly give the parity, ((—1)™), which may 
also be obtained from the P m [106-110]. Most informa- 
tively, one can use RF spectroscopy to directly measure 
P m [134], in particular by quenching the state of interest 
to a deep lattice. P m could also be measured using in- 
teraction blockade in an optical superlattice [135]. More- 
over, it may be possible to extend in- situ single-site res- 
olution experiments capabilities to measure P m directly 
[136]. 

In Fig. 8 we show our DMRG results for systems with 
open boundary conditions and L < 24 for N < 5. More 
specifically, we present the average over all sites of the 
double occupancy (D), the triple occupancy (T), the 
quadruple occupancy (Qa) and of the quintuple occu- 
pancy (Q5), defined in Appendix B. As can be seen, the 
results do not depend strongly on N. In particular in 
the limit of large U/t, we observe that the results for 
N = 4 and N — 5 are practically indistinguishable. In 
this limit, the quantities follow a power law ~ (U/t) 71 , 
with exponents (obtained for N = 5) rjn —1.9, rjx ~ 
-4.1, VQi » -7.9, 7 lQs a -12.0 at filling n = 1/(2JV) 
and r] D « -2.0, r, T « -4.0, i]q 4 « -7.9, 7/q 5 « -11.6 
at filling n = l/N. Since the double occupancy is the 
largest quantity, we analyze it in more detail. The in- 
sets of Fig. 8 show that the behavior at fillings 1/(2N) 
and l/N differs at small values of U/t. While (D)(U) 
appears to decay monotonically for all values of U/t at 
filling l/(2N), at filling l/N we identify for the small sys- 
tems a rounding or a peak for N > 2. This might be due 
to the metal-insulator transition. However, the positions 
of these maxima do not coincide with the minima of the 
fidelity susceptibility, so that we conclude that addressing 
the double occupancy in experiments with small systems 
is not sufficient to locate the phase transition. This can 
be understood since for the BKT transition all deriva- 
tives of the energy as a function of U/t behave regularly. 
According to the Hellmann-Feynman theorem, the dou- 
ble occupancy is (to a good approximation for N > 2) the 
first derivative of the energy with U, so that it is not ex- 
pected to show singular behavior at the Mott-transition 
in these systems. 

However, all these quantities show a crossover to the 
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FIG. 8. (Color online) Site averaged double, triple, quadru- 
ple and quintuple occupancies as a function of U/t for small 
systems at filling (a) n = 1/(2N), and (b) n = l/N. The 
system sizes in (a) are L — 24 for N = 2, 3, 4, and L = 20 
for N = 5. The system sizes in (b) are L = 20 for N = 2, 4, 
L = 21 for N = 3, and L = 10 for N = 5. The black straight 
lines are guide to the eyes showing the power-law behavior. 
The insets show the double occupancy as a function of U/t 
for the different values of N. 

aforementioned power-law behavior at values of U ~ lOt. 
Toghether with the behavior of the energy, this further 
supports that the minimal value of U/t for emulating 
Heisenberg physics to < 1% accuracy is approximately 
10. We will now turn to the behavior of the correlation 
functions which further support this result. 



C. Correlation functions 

1. Bosonization results for the correlation functions and 
Luttinger parameters 

According to Ref. 5, and in agreement with our results 
for the fidelity susceptibility discussed in Sec. IV, the sys- 
tem for N > 2 is in a metallic (LL) phase at small, but 
finite values of U/t. Bosonization shows that at low en- 
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ergy the spin and charge degrees of freedom separate, and 
both sectors are described by the Luttinger liquids with 
corresponding Luttinger parameters K p (charge) and K a 
(spin) [137]. Here, K a = 1 for any value of U. The lead- 
ing order contributions to various correlation functions 
can be obtained from standard Abelian bosonization [5]. 
We obtain for the density-density correlations [defined in 
Eq. (10)] 



(jV total ( r )7V total (0)) 



NK n 



2 (Trr) 



cos (2kpr) 



- r 2K p /N+2-2/N ' 

(21) 

where kp = nir/N, and n the density. For the spin-spin 
correlations [Eqs. (7) and (6)] we obtain for a ^ (3 



<S£(r)S£(0)) = - 



and for a — /3, 



1 



2(vrr) 



cos (2kpr) 

r 2K„/N+2-2/N • 



(22) 



(SS(r)SS(0)) 



N 



(K p — 1) /N + 1 B[ cos (2k F r) 



2 (Trr)" 



1 r 2K p /N+2-2/N ' 

(23) 

Note that due to the SU(iV) symmetry, Eq. (22) is, up to 
a factor of 2, the same as (S z (r)S z (0)). Also note that 
we neglect possible Akp contributions and logarithmic 
corrections in the above expressions. At U = 0, K p = 1, 
and as the repulsive interaction U is increased, the charge 
Luttinger parameter gradually decreases. As discussed in 
Ref. 5, at a sufficiently large value of U = U c , the multi- 
particle umklapp scattering terms become relevant and a 
charge gap opens, leading to the metal-insulator transi- 
tion of the BKT type discussed in the previous sections. 
In the Mott-insulating phase U > U c , the spin correla- 
tions are then simply obtained from Eqs. (22) and (23) 
by setting K p = 0, and the density correlations decay 
exponentially. 

Expression (21) can be used to obtain K p numerically. 
In the limit k — > the charge structure factor behaves as 



Af{k -> 0) 



NK, 



2?r 



£ |fc|; 



(24) 



Kp consequently can be determined by fitting the slope 
of the numerically obtained N(k) in the vicinity of k = 0. 



2. DMRG results 

In Fig. 9 we compare DMRG results for the spin cor- 
relation functions [Eq. (7)] of N = 2,3, and 4 SU(iV) 
Hubbard chains at unit filling for U — 2t, 8t, 12t to the 
spin correlation functions of the corresponding Heisen- 
berg chains [Eq. (6)]. Already at U = 2t the Heisenberg 
model reproduces the qualitative features (algebraic de- 
cay and 2k p oscillations ) of the Hubbard model. How- 
ever, the difference in the actual values shows that this 
value of U /t is outside the quantitative regime of validity 
of the Heisenberg model. For U = 8t and 12t, however, 
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FIG. 9. (Color online) Comparison of spin correlation func- 
tions [Eqs. (7) and (6)] of the 1/N filled SU(iV) Hubbard and 
SU(iV) Heisenberg chains at U = 2t, St, 12t for N = 2, 3, 4. 



the agreement is quantitative for the three values of N 
shown. The largest difference is in the nearest-neighbor 
correlations. Upon increasing N the difference decreases, 
corroborated by computing the distance between the spin 
correlation function of the Hubbard systems [Eq. (7)] and 
of the Heisenberg systems [Eq. (6)] which we define as 



d = JE [S(r)-S"(r)f 



(25) 



In Fig. 10, we see that this distance decreases with 
increasing U/t and N. For U > I2t, we find d < 0.01 for 
all values of N. Note that this criterion is matched for 
smaller values of U jt when increasing N. 
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FIG. 10. (Color online) Distance d [Eq. (25)] between the 
spin correlation functions of the l/N filled Hubbard and of 
the Heisenberg chains [Eqs. (7) and (6)] as a function of U/t 
for TV < 4. 



Figure 11 presents results for the spin correlation func- 
tions of the SU(iV) Hubbard chains at 1/2N filling. The 
results show the same characteristics as at 1 /N filling. In- 
terestingly, although this cannot be mapped to a Heisen- 
berg model, the results for U = 8t and U — 12t are very 
similar to each other for the displayed values of N. This 
suggests that the behavior in this region might be gov- 
erned by SU(iV) t— J models. These effective models cap- 
ture the interplay of spin-exchange interactions expected 
for large values of U/t with the electron itineracy. In the 
SU(2) case, it is known that this model possesses a rich 
phase diagram with superconducting phases [138]. In the 
SU(iV) case, the question arises if the phase diagrams of 
these models at N = 2 and N > 2 remain similar, as in 
the case of unit filled Hubbard chains, or if the enhanced 
symmetry might lead to unconventional phases, e.g., ex- 
otic singlet-superconductivity with singlets formed by N 
particles. 



D. Structure factors and Luttinger parameter 

Figure 12 shows our results for the various structure 
factors defined in Sec. HE at U — t and U = 15t. As 
expected from bosonization (see Sec. VC 1), all structure 
factors show a peak or a shoulder at 2/cf originating from 
the oscillatory component of the correlation functions. 
At U = t, the behavior of all structure factors at small 
k is linear with k up to k ~ 2kp (only Af a ,p(k) shows a 
nonlinear behavior). The momentum distribution func- 
tion indicates the presence of a discontinuity; this is an 
artifact due to the small system sizes available, and for 
N > 2 one would obtain a singularity in the derivative of 
n(k) at kF> according to LL theory. For N = 2, the re- 
sults look similar due to the pronounced finite size effects 
caused by the exponentially slow opening of the charge 
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FIG. 11. (Color online) Spin correlation functions [Eq. (7)] 
of SU(JV) Hubbard chains at half filling n — 0.5 for the same 
parameters as in Fig. 9. 



gap. This is also the reason why Af(k) appears to be 
linear for N = 2 despite the presence of a charge gap. 
Below we will exploit the fact that for U / 1 small enough, 
N(k) behaves linearly and obtain K p {U) from Eq. (24). 

For U = 15t, deep in the Mott-insulating phase, the 
linear behavior at small k is less pronounced or absent 
due to the finite charge gap, which leads to an exponen- 
tial decay of the correlation functions. The most drastic 
changes are seen in Af(k) and n(k), which directly probe 
charge degrees of freedom. In these quantities, the sin- 
gularities at 2kF and kp disappear, as expected for a 
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Mott insulator with a large charge gap. Note that S(k) 
behaves linearly in the region < k < 7r/4, and is the 
same for all values of N with a slope of l/(27r), in agree- 
ment with the bosonization result for the spin correlation 
function, Eq. (22). 

Fig. 13 shows our results for K p {U) obtained from fit- 
ting the slope of N~(k) at k « 0. Note that for U > U c , 
formally K p does not enter the correlation functions [137] 
and the structure factor cannot be used to determine K p . 
However, due to the exponentially slow opening of the 
gap, when the system size is much smaller than the cor- 
relation length the structure factor for the finite systems 
appears to behave linearly so that we fit the slope also in 
these cases. It appears that for N = 2, there seems to be 
an inflection point at U = 0, which leads to a minimum of 
x(U) computed from K p using Eq. (15). For N > 2, sim- 
ilar inflection points seem to appear at U ~ 1.5t(N=3) 
and U ~ 2i(N=4), i.e., close to the values of J7 m i n at 
which the fidelity susceptibility x(J7) has its minimum. 
Additional inflection points seem to appear at larger val- 
ues of U (N = 2: U w 2i; N = 3: U « At; N = A: 
U ~ At). These are in rough agreement with the position 
of the maxima of x(U) discussed in Sec. IV C. 

VI. SUMMARY AND CONCLUSIONS 

Using BA, bosonization, and DMRG, we have inves- 
tigated SU(iV) Hubbard chains and identified the region 
of validity of both the strong- and weak-coupling pertur- 
bative regimes. For N > 2, where the BA is known to be 
an approximation, the values of the energies for N = A 
agree with the ones obtained by DMRG with a relative 
error < 4%. We therefore use the BA to explore the be- 
havior at large N which is difficult to access with DMRG. 
In addition, by computing the fidelity susceptibility, we 



have shed new light on the value of the critical interac- 
tion strength U c for the Mott transition at 1/N filling. 
We identify a clear minimum in the fidelity susceptibility 
in the vicinity of the putative U c even for rather small 
system sizes. Since the same behavior is obtained for 
N = 2, for which U c = is known exactly, we conclude 
that U c > t for all N > 2. For experiments with N > 2, 
this signifies that it should be possible to observe Lut- 
tinger liquid behavior for U < t even at 1/N filling. 

For large U/t, we identify that SU(iV) Heisenberg mod- 
els provide a very accurate (< 1% error) description of 
the system for U > 12t for all values of N. For these 
values of U /t, the absolute difference between the corre- 
lation functions of the Hubbard and the Heisenberg sys- 
tems is < 0.01. We expect therefore that in this regime 
the SU(iV) Heisenberg model will quantitatively describe 
experiments with ultracold alkaline earth atoms on op- 
tical lattices. We expect that also in higher dimensions 
this may be true for similar values of U/t. Given that it 
is more favorable for the experiments to work with values 
of the spin-exchange interaction J as large as possible, we 
therefore suggest to search for the proposed chiral spin 
liquid state [67] in experiments on square lattices for val- 
ues of the interaction lOt < U < 15t. 
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Appendix A: Bethe ansatz equations 



Here we summarize the basic approximate Bethe ansatz equations developed in Refs. 94 and 95. 



The rapidity distributions p c are associated with the charge degree of freedom, and the pj's for j e {1, . . . , s} are 
associated with the rapidity distribution governing the difference in spin states j — 1 and j (where we interpret spin 
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U =15 




FIG. 12. (Color online) Structure factors of the various correlation functions Eqs. (7) - (11). (a) Charge structure factor 
[Eq. (10)] for U = t. (b) Spin structure factor [Eq. (7)] for U = t. (c) Structure factor Af a , a (k) [Eq. (8)] for U = t. (d) 
Structure factor J\f a ,p{k) [Eq. (9)] for U = t. (e) Momentum distribution function [Eq. (11)] for U — t. (f)-(j): the same 
quantities, but for U — 15i. 
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state "0" as the charge). These distributions satisfy the coupled set of N linear integral equations 



/'-</•■> = ^ + "~7r~~^ dAK^sink- A) Pl (A) 
_'m Ztt ./_ Ai 



PiW = 7^l dk K^A- sin k)p c (k)- — I dA' K 2 {A-A')p x {A') 

dA' K 1 (A-A') P2 {A') 

2 

P,[A ) = - / <LV K,(A - A')pa-i(A') - ~ y dA'K 2 (A - A>,(A') 

dA' X! (A - A')p s+ i (A') for s = 2, . . . , N - 2, 



2tt 7 j. v /r v ' 2tt 

1 



1 



2tt 



-A, +1 



i /"Ajv-2 i /*Ajv— i 

p N . 1 (A) = ~ dMKx^k-M^n-^M)-— I dA'K 2 (A-A') PN ^(A') (Al) 

2tt '-a n _ 2 2tt J^^j 



with 



The parameters fco and A s for s = 1, . . . , N — 1 are determined by the charge and spin densities n c and nj through 

dk p c (k) 

— ko 

/ho i'Ax 
dkp c (k) - I dA Pl {A) 
-ko J —Ax 

dA/3 a _i(A) 

A a -i 

dAp s (A) for s = 2, . .. , AT - 1 (A3) 

'-A, 

and the ground state energy per site is given by 

dk cos(k)p c {k) (A4) 

-ko 

We numerically solve the integral equations, Gauss-Legendre discretizing the linear integral equations and solving 
the resulting linear equations [139]. To apply the discretization procedure for finite intervals, we first transform 
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the intervals (— fco>^o) an d (— A S ,A S ) to ( — 1,1). A uniform rescaling of the coordinates is suitable for the (— ko,ko) 
interval, but is unfavorable for the rest because A s tends to infinity for the population balanced gas, while the width of 
the function in the original units approaches a constant. Thus, a uniform rescaling would require an unnecessarily large 
number of points as one would sample mostly where the integrand was zero. Instead, we rescale to new coordinates 
defined by 

"(A) = (A5) 

Defining the Jacobian 

• i \ a 1 + A + Au 2 

JA(U) = A (1 + A-A U ^' ( A6 ) 



the Bcthe ansatz equations Eq. (Al) become 

. 1 cos(fc o) Z" 1 , . . . T . ( , , a(Ax)u \ 
*(«) = ^ + J_duj Al (u) Kl (sMkoq) - 1 a{Ai)u2 ) *(«) 

*(«) = | / * *i (i^A^a - sin ( fc o?)) *(?) 

2tt \ v 1 ^ a ( A i) u 1 - a(Ai)(u') 2 / 

1 /" , / • ( i\ v ( a ( A i) u a{A 2 )u' \ 

27r j \ v l-a(Ai)u i; 1 - a{K 2 ){u'Y J 

27T J_x \l-a(A s )u 2 1 - a(A s _i)(u') 2 / 

- 7T / du 3hK 2 r 2 - ,. 2 Psiu ) 

27r y_! v 1 ^"!^) - " i - Qt(Aa)(u') / 



a{A s )u a(A 8+ i)u' 



2tt y_ x J ' ls+1 'Vl-QfA^ l-a(A s+1 )(u') 2 



Pa+i(«') for s = 2,...,iV-2, 



PN-l(U) = — du JK N _ 2 K X - — r-g- - — -r— ^ M ) 

2ttJ_ 1 VI - a(A W -i)u^ 1 - a(A A r_ 2 )(w ) 2 



27r y_x \1 - a{A N -i)u- 1 - a(A N -i)(u'y J 



2tt 

The relevant charge density, spin densities, and ground state energy are given by 



n c = k J dqp c (q) 

m = k J dqp c (q)- J duj Al (u)p 1 (u) 
n s = J dujA^iujps-xiu) 

du jA s (u)p s (u) for s = 2, . . . , N — 1, 

-1 



£ba = -2tk J dqp c (q). 



(A8) 



Note that the BA results are obtained in the thermodynamic limit. In Sec. Ill, we compare the results obtained by 
this procedure to the umerical results of the DMRG after extrapolating to the thermodynamic limit. 
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TABLE II. Numerical values of the coefficients in Eqs. (B1)-(B4). 



Appendix B: Expressions for the higher local occupancies 



In this appendix we define the expressions for the double occupancy (£>), the triple occupancy (T), the quadruple 
occupancy (Q4), and the quintuple occupancy (Q5) discussed in Sec. VB. We obtain: 



(Qa 



(Qs) 



1 L 

i=i 
1 L 



-If Y, («'"?"> -rf Y, (nfnf'<"<"')-r 5 D «n?n?nJn|) 



a. a ,a ,a 



i=l 
L 



a, a .a .a 



i=l 



£ («'<"<"') - if* (n* n? n! n\ nf ) 



a, a' .a" ,a'' 



(Bl) 

(B2) 

(B3) 
(B4) 



The sums over the flavors a, are over all possible permutations, and the numerical coefficients T° are listed in Tab. II. 
To obtain these coefficients, we choose them so that the N'th order polynomial of nf reproduces the action of P m on 
a complete basis: f(n a ) \m') = P m \m') — S mm i for the N + 1 values of m = 0, . . . , N. 
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